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Abstract 

The computer code DECAY4 is developed to generate initial energy, time and angular 

distributions of particles emitted in radioactive decays of nuclides and nuclear (atomic) 

deexcitations. Data for description of nuclear and atomic decay schemes are taken from 

the ENSDF and EADL database libraries. The examples of use of the DECAY4 code in 

O ' several underground experiments are described. 
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^\ 1 Introduction and overall description of the DECAY4 
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Despite the fact that effect to background ratio is the central problem of all experimental 
physics, there is a certain class of experiments for which this problem is so crucial that even 
the possibility of their performance itself depends strongly on the reached background level of 
the used detectors. These are so called underground experiments devoted to investigation of 
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the very rare or forbidden decays and processes like, for instance, double beta decay, proton 
decay, dark matter particles search, solar neutrino study and so on. The ultimate sensitivity of 
such experiments is determined mainly (except the available source strengths) by the detector 
background. The first origin of background is due to cosmic rays and can be eliminated by the 
proper underground site for the set up. The second (and the most crucial for sensitivity) source 
D . of the background is the decays of the nuclides from the radioactive impurities in the detector 
itself, in the materials used for detector mounting and shielding, and in the surroundings. 
Therefore it is obvious that simulation of the background and, in particular, simulation of 
the nuclides decays is the overwhelmingly important part of such kind of research which can 
allow: i) to understand and determine the origins of the background and hence to find certain 
methods to eliminate or suppress the background contributions; ii) to build up the background 
model and response functions of the detector for the effect being sought (together with the 
detector energy and efficiency calibrations, resolution, source activities, etc.), thus to extract 
and evaluate searched effect (or to exclude it) more precisely. 

There are several general programs, which are commonly used for simulation of the particles 
interactions in the experimental set up as, for example, GEANT package |I| or EGS4 code 0. In 
any such a program user should describe initial kinematics of events by using the so called event 
generator. The last is an important part of the simulation program providing the information, 
which particles and how many of them are emitted, what are their energies, directions of 
movement and times of emission. Existing computer codes RADLST || and IMRDEC [[| only 
determine the radiation spectra due to decay of nuclides, and therefore could not be used for 
the further particles tracking. 

In attempt to cover this lack, the code, named DECAY4, was developed for generation of 
events in the low energy nuclear and particle physics (double beta decay and decay of radioactive 
nuclides). This code was elaborated during the last decade, mainly for 2(3 decay research [[J. 
The plan of the paper is as following: first, the overall features of the DECAY4 generator and 
used databases are considered, then detail descriptions of the parts associated with the double 
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beta decay of atomic nuclei and decays of the natural and artificial radioactive nuclides are 
presented. In the last section several examples illustrating the use of the DECAY4 generator 
in the real underground experiments are shown. 

The developed program DECAY4 gives the possibility to generate events of the 2(3 decay of 
atomic nuclei and of the radioactive decays (a, P ± , p, n decay, electron capture) of all known 
unstable isotopes. It is divided into two main sections: 

a) INIT - search and reading of all parameters of the nucleus and its decay needed for decay 
simulation from the ENSDF || (or NuDat 0), EADL || and other libraries, in order to build 
up the nuclear and atomic decay schemes; 

b) GENDEC - a Monte Carlo events generator itself. 

The ENSDF database library includes the following information on about 2500 isotopes 
used for generation of the radioactive decays: a) decay modes, their probabilities and energy 
releases, isotopes half-lives; b) radiation type, particles energies and intensities; c) parameters of 
nuclear levels (half-life, spin, parity and excitation energy); d) parameters of nuclear transitions 
(branching ratios, multipolarities, coefficients of internal conversion and mixing ratios). 

The DECAY4 code uses also tables with the data on atomic properties of isotopes (electron 
binding energies, electron capture (EC) subshell ratios, X rays and Auger electrons intensities) 
from the EADL database || (as well as from |J) and tables with theoretical Hager-Seltzer 



conversion coefficients 10 



The GENDEC part of the DECAY4 generates the energy, time of emission, direction and 
polarization for the following emitted particles: 1) electrons and positrons from single and 
double (3 decay; 2) a particles from a decay, protons and neutrons from p and n decays; 3) 
7 quanta from nuclear deexcitation process; 4) conversion electrons; 5) e~e + pairs from the 
internal pair conversion; 6) 7 quanta from bremsstrahlung in (3 decay and EC; 7) neutrinos 
(antineutrinos) from EC or (3 (2/3) decay; 8) X rays and Auger electrons from the atomic 
deexcitation process. 



2 Double beta decay processes 

The DECAY4 describes double beta processes (2(3" and 2(3 + decays, electron capture with 
emission of positron e(3 + and double electron capture 2s) for all nuclides. 2(3 transitions to the 
ground state as well as to excited + and 2 + levels of the daughter nucleus are allowed. If 2(3 
decay occurs to an excited level of a nucleus, the electromagnetic deexcitation process follows. 
The energy release Qpp for double beta processes is taken from the table of the atomic masses 
nj . For each transition to the ground or an excited level, various modes (with emission of 
two neutrinos or Majoron, neutrinoless decays due to nonzero neutrino mass or right-handed 
admixture in the weak interaction, etc.) and mechanisms (two-nucleon 2n and A-isobar N*) of 
double beta decay are possible. Below we give the list of double beta processes which could be 
simulated with the DECAY4: 

1. 0z/2/5 ± decay with neutrino mass, + — + transition, 2n- mechanism; 

2. §v2(3 ± decay with right-handed currents, + — + transition, 2n-mechanism; 

3. Qv2(3^ decay with right-handed currents, + — + and + — 2 + transitions, N*-mechanism; 

4. 2u2(3 ± decay, + — + transition, 2n-mechanism; 

5. 0^2/3 ± decay with emission of Majoron, + — + transition, 2n-mechanism; 
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6. 0z/2 / 9 ± decay with double Majoron emission, + — + transition, 2n-mechanism; decay 
with charged L = —2 Majoron or massive vector Majoron; 

7. 0z/2/5 ± decay with right-handed currents, + — 2 + transition, 2n-mechanism; 

8. 2z/2/3 ± decay, + — 2 + transition, 2n- and N*-mechanisms; 

9. 0ve/3 + decay; 

10. 2ve(3 + decay, + — + and + — 2 + transitions; 

11. 0v2e decay; 

12. 2v2e decay, 0+ - 0+ and 0+ - 2+ transitions. 

The theoretical formulae for the energy and angular distribution p(E\, E 2 , cosO) of emitted 
electrons or positrons are based on works |12|, |H| [14], [15|, For example, for the first process: 

p(E x , E 2 , cos9) = p l {E 1 + 1)F(E U Z) p 2 {E 2 + 1)F(E 2 , Z)S(E -E x - E 2 ){1 - fcfacoaB), (1) 

where Ei is the kinetic energy of the z-th e~ or e + (in units of the electron mass m e c 2 ), pi is 
the momentum (in units of m e c), F(Ei,Z) is the Fermi function, Z is the atomic number of 
the daughter nucleus (Z > for 2(3~ and Z < for 2(3 + decay), 6 is the angle between the 
particles directions, E is the energy available for the particles (E = Qpp — E? x for 2(3~ decay 
and E = Qpp — 4 — E ex for 2f3 + decay, E ex is the energy of the populated level of the daughter 
nucleus), $ = Pi/(Ei + 1). 



3 Radioactive decays of nuclides 

The DECAY4 describes six decay modes: (3~, a, p and n decays, electron capture and 
f3 + decay (EC) and isomeric transition (IT). The modes d (d = /3~ , a, p, n, EC, IT), their 
probabilities p d , available decay energies Q d and isotopes half-life T\/ 2 are taken from the 
ENSDF or NuDat databases. The decay mode d is sampled according to the probabilities 
p d . 

(3~ decay. The endpoint energies in P~ decay E° are related with the energy release 
and the level energies E\ x of the daughter nucleus by the equation 

QP- = E° + E\ x . (2) 

Kinetic energy of the beta particle E is sampled in accordance with the distribution 

p{E)=p.(E + l). (£° - E) 2 ■ F(E, Z) ■ S k (E) , (3) 

where Sk {E) is the factor of forbiddenness. Probability of the internal bremsstrahlung in beta 
decay and energy-angular distribution of bremsstrahlung 7 quanta are calculated as in ]T7[ . 

a, p, n decays. The energies of particles E\ (k = a, p, n) are related with the level energies 
Ef x of the daughter nucleus by the equation 

E* = A d /A p -(Q k -En, (4) 

where A p , is the mass numbers of the parent and daughter nuclei, respectively. 

EC (electron capture and (3 + decay). The ENSDF database includes the information 
on the probabilities pf c (for EC) and pf (for /3 + decay) for i-th level of the daughter nucleus. 
If the level is populated in j3 + decay, energy of positron is sampled according to (3), where 

£° = Q EC - 2 - Ef. (5) 



3 



If the i-th nuclear level was populated in the electron capture, the number of the atomic 
subshell x (x — K, Li, L 2j L 3; Mi, M 2 , M 5 , Ni, N 2 , N 7 , Oi, 2 , 7 ), where primary 
electron vacancy is created, is sampled according to probabilities P x c |§ : 

P X EC (Z, q x ) = const ■ n x p 2 ^ ^q 2 J L - k ^ (3 2 x B x /[(2k x - 1)!(2L - 2k x + 1)!], (6) 

where L is the electron capture transition angular momentum, n x is the relative occupation 
number for partially filled subshells x (n x = A^/iV™ 3 *, N x is the number of electrons in the 
subshell x, iV™ ax is the maximal number of electrons in the subshell), q x = Q EC — E^ x — E x 
is the neutrino energy, E x is the electron binding energy in the parent atom, and k x is x 
subshell angular momentum. The squared amplitudes 0l,B x p^ x ~^ of the bound-state electron 



radial wave functions and the electron binding energy E x are taken from |18], ||. The internal 
bremsstrahlung probability and spectra of bremsstrahlung 7 quanta in allowed electron capture 
transition from x atomic subshell is calculated in accordance with 



Nuclear deexcitation process. The nuclear deexcitation process occurs if a daughter 
nucleus is in the excited z-th level with energy Ef x . The electromagnetic transition from i-th 
to j-th level is sampled according to probabilities 

Pij = Jijl "A?' (^) 

where Jij is the branching ratio of electromagnetic transition from i-th level to j-th taken from 
the NuDat or ENSDF databases. 

There are three possible modes of electromagnetic transition with emission of: (1)7 quantum 
with energy E 1 = Ef x — E ex \ (2) conversion electron with energy E x e = E^ x — E ex — E x (the 
condition E x e > should be fulfilled), here E x is the electron binding energy on the x subshell; 
(3) conversion electron-positron pair with total energy E cp = E\ x — E ex — 2 (if E cp > 0). To 
sample the mode, the respective probabilities p]j, p^, are used, where 

pi = 1/(1 + afj + of), p% = a%pl, p% = a%pl (8) 



for all transitions beside E0 



pi = 0, pt 5 = 1/(1 + p% = ■ p% (9) 

t ce a cp 



for E0 transition, where af? 0$ are the coefficients of the internal electron and pair conversion, 



respectively, and I? are the intensities of internal electron and pair conversion, respectively. 

Total a^, partial subshell a^(s m ), shell a^(s) coefficients of electron internal conversion is 
related by 

< = E s W; «S00 = E m (10) 

where s is the shell index and m is the subshell index of the s shell (s m = x). Coefficients 
and a^ e (s) are taken from the NuDat or ENSDF databases. In case if their values are absent, 
the coefficients are calculated as 

a%{s m ) = [a c l ](s m ,E\n 1 Xi,Z) + 5 lj a^(s m ,E\n 2 X 2 ,Z)]/(l + S 2 J ), (11) 

where the values of partial coefficients a ce (s m , E 1 , tt\, Z) are taken from [jnj. Here ir\ is 
transition multipolarity, and 5ij is the mixing ratio of different multipolarities in the i —> j 
transition. For pure E0 transitions, electron conversion coefficients are calculated according to 
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formulae from |J. Coefficient of internal pair conversion a c ? is given by formulae similar to (10), 
and partial coefficients a cp (i? 7 , 7rA, Z), and electron-positron energy and angular distributions 
are calculated according to formulae from ||, [19 . 



X rays and Auger electrons. The vacancies in the atomic shells are created in the 
electron capture or in the result of internal electron conversion. Ionized atom is deexcited 
filling the vacancies by electrons from higher atomic shells; X rays and Auger electrons are 
emitted in this process. For the sampling of the atomic deexcitation process, the values of the 
electron binding energies and occupation numbers, radiative and radiationless partial widths are 
taken from the EADL library ||. The type of process (radiation or Auger electron emission) 
is sampled accordingly to the values of radiative and radiationless partial widths. For the 
vacancy in the Sj subshell, the energy of X ray for radiative process, in which higher vacancy 
qj is created, and the energy of Auger electron for nonradiative process (vacancies qj and U are 
created) are given by 

E*_ q . = E St (Z) - E q .(Z); E*_ q . ti = E St (Z) - E Qj (Z) - E tl (Z) - AE q , tr (12) 

Correction AE qjtl is found from equations in pD| . Multivacancy corrections for energies and 
partial widths are included also. The X ray and Auger electron energy spreading is accounted 



according to Lorentzian distribution function |20 



Decay chains and time characteristics. Time interval T k between appearance of the 
unstable state k of the nucleus or atom and its decay with emission of the particle is sampled 
in accordance with the relation 

T k = -T} /2 ]n V /]n2, (13) 

where 77 is a random number uniformly distributed in the range (0,1), Ty 2 is a half-life of 
the unstable state. The DECAY4 can also simulate the full decay chain of the parent nuclide 
together with its daughters. As an example, the generated spectra of particles, emitted in 
decay of nuclides from full 238 U chain, are shown in fig. 1. Activities of all daughter nuclei are 
calculated by the DECAY4. 

Angular correlations between emitted particles. The direction (polar and azimutal 
angles 9i, (pi) of the particle i is sampled as isotropical if one of the next conditions is fulfilled: 
1) the particle is the first in the decay of unoriented parent nuclear state i; 2) ij < 1 (ij is 
nuclear spin before emission of the particle i); 3) T % > r max (T l is the time given by (13), r max 
is user defined parameter (time) to account the influence of external fields which violate the 
angular correlation); 4) \ = or A,_i = (Aj and A,_i are angular momenta of particles i and 
i — 1, if particle i — 1 exists); 5) one of the values Ii-i, ij, ij+i, Aj, Aj_i is unknown and could 
not be evaluated. In case, if particle % + 1 does not satisfy any of the above conditions but 
particle i satisfies one of them, the polar angle 9\ +l of particle % + 1 in the coordinate system 
of particle % is sampled in according to correlation function [21] 



u 

W(^ = £^; =0 ^ i+1 (AiA^^ (14) 

where A- = \ ± 1 is the second possible angular momentum of the particle (for the mixture of 
multipolarities), Xi is type of the particle i (xi = 7, f3, a, e), v4^ +i (AjA-/j + i/j, Xj) are angular 
momenta functions, and k is even. In this case correlation function is independent from the 
azimutal angle <p\ + i- If the next emitted particle i+2 exists, which does not satisfy any condition 
among l)-5), its direction is determined by the more complicated correlation function, which 
depends on directions of previous particles. Correlation of linear polarisations of 7 quanta is 



max 



also taken into account [21 
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Figure 1: Generated initial energy spectra of particles emitted in the decay of 238 U chain (in 
equilibrium): (a) 7 quanta, (b) electrons, (c) a particles, (d) antineutrinos. 

4 Conclusions 

The code DECAY4 was successfully used in several underground experiments for detectors 
design and optimization, simulation of the backgrounds (with the help of the GEANT package) 
and for evaluation of the results. Some examples are listed below. 

1. Kiev 2(3 decay experiments performed in the Solotvina Underground Laboratory in a salt 
mine 430 m underground (~1000 m w. e.): 

1.1. Cadmium tungstate crystal scintillators (enriched in 116 Cd to 83%) were used to study 
116 Cd p2| , p3| . The background of the 116 CdW04 crystal (15.2 cm 3 ) in the energy region of 
interest (Q / 3 j a=2805 keV) was equal to ^0.6 counts/y-kg-keV. With the 19175 h statistics the 
half-life limit for neutrinoless 2(3 decay of 116 Cd was set: T / /2 > 3.2 • 10 22 y (90% C.L.) [f£J. 
Limits on Oz/ modes with emission of one (Ml) or two (M2) Majorons were established also: 
T^f 1 > 1.2 • 10 21 y and T^f 2 > 2.6 • 10 20 y (90% C.L.) 0. Comparing these limits with the 
theory, the restrictions on the neutrino mass (m u ) < 3.9 eV and Majoron- neutrino coupling 
constant g < 2.1 • 10~ 4 were derived, which are among the most sensitive results for other nuclei 

g. 

1.2. The study of the 2/3 decay of 160 Gd was carried out with the help of the 95 cm 3 
Gd 2 Si05:Ce crystal scintillator. The background was reduced to pal.O cpd/keV-kg in the vicin- 
ity of energy (^1.73 MeV). The improved half-life limit was set for 0^2/3 decay of 160 Gd: 
T° v /2 > 1.2-10 21 y at 68% C.L. [g§ . 

2. DAMA collaboration experiments performed deep underground in the Gran Sasso Na- 
tional Laboratory: 

2.1. Two radiopure CaF 2 :Eu crystal scintillators (370 g each) have been used to study 2(3 
decay of 46 Ca and double electron capture of 40 Ca as well as for dark matter search. The 
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highest up-to-date half-life limits were reached for Oz/ and 2v double electron capture of 40 Ca: 
T^ 2 > 4.9-10 21 y and T^ 2 > 9.9-10 21 y (68% C.L.) 

2.2. The study of the 2(3 + decay of 106 Cd was performed with the help of two low background 
Nal(Tl) crystals and enriched (to 68%) 106 Cd samples (wl54 g). New T\/2 limits for the (3 + (3 + , 
/3+/EC and EC/EC decay of 106 Cd have been set in the range of (0.3-4)- 10 20 y at 90% C.L. 



3. NEMO collaboration experiment performed in the Frejus Underground Laboratory to 
study 2/3 decay of 100 Mo [E3]. For background simulation the first version of the DECAY4 code 



l29j was used. The clear two neutrino 2(3 signal (1433 events during 6140 h) was observed, 
leading to a half-life T 1/2 = 0.95 ± 0.04(stat)±0.09(syst)-10 19 y f2§. 

4. Heidelberg-Kiev collaboration: The background simulations were performed for the 
GENIUS project which aim is to increase the present sensitivity for 2(3 decay and dark matter 
search. Contributions from the cosmogenic activity produced in the Ge detectors and from their 
radioactive impurities as well as from contamination of the liquid nitrogen and other materials 
were calculated. External 7, fi and neutron backgrounds were considered also. The results of 



calculations evidently show feasibility of the GENIUS experiment ]30 . 
In result, we can conclude: 

1. The DECAY4 code is the powerful tool for simulation of the radionuclides decays in the 
wide range of decay's modes, emitting particles, etc. 

2. The DECAY4 can be connected easily with the codes which simulate the particles 
propagation, like GEANT, EGS4, etc. 

3. The code includes the most advanced databases: ENSDF, NuDat, EADL and others. 

4. The event generator DECAY4 was used successfully in many underground experiments 
(Kiev, Roma-Kiev, DAMA-, NEMO-, GENIUS-collaborations). 
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